#DESCRIPTION
#module generating slip surfaces by the slope search method as employed in Roscience' Slide

#Reference:
#Slide online manual. http://www.rocscience.com/help/slide/webhelp/Slide.htm, accessed June 18 2014

#libraries and modules
import math
import numpy as np
import matplotlib.pyplot as plt



#math/trigo functions
tan=math.tan
sin=math.sin
atan=math.atan
rad=math.radians

def slope_pts(H,B,x):
    #generating (y) points for simple slope defined by slope height (H) and gradient (B), and set of (x) points
    x=np.asarray(x)
    zero=np.zeros(len(x))
    crest=H*np.ones(len(x))
    y=np.where(x>0,x*tan(rad(B)),zero)
    y=np.where(x>H/tan(rad(B)),crest,y)
    return y

def slip_surface_generation(H,B, num_slip, min_x1_c,max_x1_c,min_x2_c,max_x2_c, min_alpha, delta_B_max_alpha,num_slices):

    #INPUTS

    #H,B:           floats, slope height (H) and gradient (B)
    #num_slip:      integer; number of slip surfaces to generate per H,B and set of material property

    #coefficients for defining endpoints relative to slope geometry:
    #min_x1_c:      float, factor of slope base length, <=0.5;                (+) within slope face; (0) toe; (-) slope base 
    #max_x1_c:      float, factor of slope base length, >=min_x1_c, <=0.5;    (+) within slope face; (0) toe; (-) slope base
    #min_x2_c:      float, factor of slope face length, >=0;                  (0) at slope crest; (+) beyond slope crest 
    #max_x2_c:      float, factor of slope face length, >=min_x2_c; 

    #min_alpha:     float; minimum exit angle of slip surface 
    #delta_B_max_alpha: float; minimum difference between B and maximum exit angle of slip surface
    #num_slices:    integer; number of slices to generate per slip surface

    #1. defining array of slip surface endpoints(x1,y1,x2,y2)
    min_x1=min_x1_c*H/tan(rad(B))
    max_x1=max_x1_c*H/tan(rad(B))
    min_x2=H/tan(rad(B)) + min_x2_c*H/sin(rad(B))
    max_x2=H/tan(rad(B)) + max_x2_c*H/sin(rad(B))

    x1=np.random.rand(num_slip*5)*(max_x1-min_x1)+min_x1   
    y1=np.where(x1<0,np.zeros(len(x1)),tan(rad(B))*x1)   
    x2=np.random.rand(num_slip*5)*(max_x2-min_x2)+min_x2
    y2=H

    #2. defining array of arc exit angles (x1alpha)
    x1alpha=np.random.rand(num_slip*5)*((B-delta_B_max_alpha)- min_alpha)+min_alpha

    #3. solving arc centers (h,k) and radii (R) from endpoints and exit angles
    m2=-(x2-x1)/(y2-y1)
    b2=(0.5*(y1+y2)) -m2*(0.5*(x1+x2))
    m1=np.tan(np.radians(270+x1alpha))
    b1=y1-m1*x1
    h=-(b1-b2)/(m1-m2)
    k=m1*h+b1
    R=np.sqrt((x1-h)**2+(y1-k)**2)

    #4. evaluating slip surfaces
    errcode=np.zeros(len(x1))
    arc_ok=0
    for i in range(len(x1)):
        if k[i]>=y2:
            xarc=np.linspace(x1[i],x2[i],5*num_slices+1)
            yarc=k[i]-np.sqrt(R[i]**2-(xarc-h[i])**2)
            yslope=slope_pts(H,B,xarc)
            if np.sum(np.where(yslope-yarc<=-0.00001,np.ones(len(xarc)),np.zeros(len(xarc))))>0:
                #ERROR slip surface intersects slope surface at points other than the arc endpoints
                errcode[i]=-105    
            else:
                #SLIP SURFACE OK
                errcode[i]=1
                arc_ok=arc_ok+1
        else:
            #ERROR reverse curvature slip surface arc center (k) < slope height (H)
            errcode[i]=-114
            
        if arc_ok>=num_slip:break
   
    return x1[np.where(errcode==1)],x2[np.where(errcode==1)], h[np.where(errcode==1)], k[np.where(errcode==1)], R[np.where(errcode==1)],errcode

def trial():
    H=50
    B=60
    num_slip=200
    min_x1_c=-1
    max_x1_c=0.5
    min_x2_c=0
    max_x2_c=0.5
    min_alpha=-30
    delta_B_max_alpha=10
    num_slices=10
    x1,x2,h,k,R,err=slip_surface_generation(H,B, num_slip, min_x1_c,max_x1_c,min_x2_c,max_x2_c, min_alpha, delta_B_max_alpha,num_slices)
    x=np.linspace(x1.min(),x2.max(),50*num_slices+1)
    y=slope_pts(50,60,x)
    
    for i in range(len(x1)):
        xarc=np.linspace(x1[i],x2[i],5*num_slices+1)
        yarc=k[i]-np.sqrt(R[i]**2-(xarc-h[i])**2)
        plt.plot(xarc,yarc,'-')
    plt.plot(x,y, 'k-',lw=3)
    plt.show()


        

    
